* Connector Paper Tables

clear all
set more off, perm
set maxvar 120000

cap log close
cd "/disk/eagedisk1/mahc_sl.work/stiebris/code_connector/"
global zip "../../Master_Formatted_Data/zip2stablewrap.dta"
global zip_income "../raw/MedianZIP-3.xlsx"
global coord "../raw/tl_2016_us_zcta510.shp"
global mail_comply "/disk/eagedisk1/mahc_sl.work/Master_Formatted_Data/JPAL_Flatfiles/Arm4_Letters/Arm4LetterSenders_unlocked.xlsx"

* Stata Packages
ssc install tabstatmat
ssc install mat2txt

program main
* LASSO Controls
	lasso_controls
* Main Tables
	balance_table
	main_results
* Appendix Table
	mm_by_zp
	het_effects
	sum_stats_mail
	rd_premium
	rd_log
end

* LASSO Controls
program lasso_controls
foreach var in "enroll90"{
use ../output/connector_regression_prep, clear

egen med_age = median(elig_age)
encode zipstring, gen(zip_factor)
gen over35 = (elig_age>35 & !mi(elig_age))
gen overmed = (elig_age > med_age & !mi(elig_age))
gen duesooner = (time_to_formdue < 30)

local select "female new_elig prior_enroll_2 spanish_letter i.WrapRegion_f time_elapsed new_fpl i.incomebracket over35 overmed wave1 duesooner i.zip3digit i.zip_factor female#new_elig female#prior_enroll_2 female#spanish_letter female#i.WrapRegion_f female#c.time_elapsed female#c.new_fpl i.female#i.incomebracket female#over35 female#overmed female#wave1 female#duesooner female#i.zip3digit female#i.zip_factor new_elig#prior_enroll_2 new_elig#spanish_letter new_elig#i.WrapRegion_f new_elig#c.time_elapsed new_elig#c.new_fpl new_elig#i.incomebracket new_elig#over35 new_elig#overmed new_elig#wave1 new_elig#duesooner new_elig#i.zip3digit new_elig#i.zip_factor prior_enroll_2#spanish_letter prior_enroll_2#i.WrapRegion_f prior_enroll_2#c.time_elapsed prior_enroll_2#c.new_fpl prior_enroll_2#i.incomebracket prior_enroll_2#over35 prior_enroll_2#overmed prior_enroll_2#wave1 prior_enroll_2#duesooner prior_enroll_2#i.zip3digit prior_enroll_2#i.zip_factor spanish_letter#i.WrapRegion_f spanish_letter#c.time_elapsed spanish_letter#c.new_fpl spanish_letter#i.incomebracket spanish_letter#over35 spanish_letter#overmed spanish_letter#wave1 spanish_letter#duesooner spanish_letter#i.zip3digit spanish_letter#i.zip_factor i.WrapRegion_f#c.time_elapsed i.WrapRegion_f#c.new_fpl i.WrapRegion_f#i.incomebracket i.WrapRegion_f#over35 i.WrapRegion_f#overmed i.WrapRegion_f#wave1 i.WrapRegion_f#duesooner i.WrapRegion_f#i.zip3digit i.WrapRegion_f#i.zip_factor c.time_elapsed#c.new_fpl c.time_elapsed#i.incomebracket c.time_elapsed#over35 c.time_elapsed#overmed c.time_elapsed#wave1 c.time_elapsed#duesooner c.time_elapsed#i.zip3digit c.time_elapsed#i.zip_factor c.new_fpl#i.incomebracket c.new_fpl#over35 c.new_fpl#overmed c.new_fpl#wave1 c.new_fpl#duesooner c.new_fpl#i.zip3digit c.new_fpl#i.zip_factor i.incomebracket#over35 i.incomebracket#overmed i.incomebracket#wave1  i.incomebracket#duesooner i.incomebracket#i.zip3digit i.incomebracket#i.zip_factor over35#overmed over35#wave1 over35#duesooner over35#i.zip3digit over35#i.zip_factor overmed#wave1 overmed#duesooner overmed#i.zip3digit overmed#i.zip_factor wave1#duesooner wave1#i.zip3digit wave1#i.zip_factor duesooner#i.zip3digit duesooner#i.zip_factor i.zip3digit#i.zip_factor"

splitsample, gen(sample) split(0.75 0.25) rseed(300)

qui: lasso linear enroll90 (arm2 arm3 arm4 i.batch) `select' if use_id_90==1 & sample==1, rseed(330)
est sto lambda
lassocoef lambda, sort(coef, standardized) nofvlabel
di "`var' lasso results"
}
end

* Main Tables
program balance_table
* Balance Table
use ../output/jpal_sample, clear
keep if use_id==1 & letter_mailed < mdy(7,1,2019)
sort arm 

* Summary Statistics by Arm
eststo clear
by arm: eststo: estpost summarize ///
    elig_age elig_fpl female prior_enroll_2 zeropremium_eligible new_fpl new_elig

esttab using "../output/main_results/tables/jpal_prelim_summarytable.tex", replace ///
cells("mean sd")  label ///
title("Descriptive statistics for full study population") ///
mtitle("Arm 1" "Arm 2" "Arm 3" "Arm 4") tex

forval i = 1/7{
	gen batch`i' = batch == `i'
	}

drop batch
forval i = 1/5{
gen incomebracket`i' = incomebracket == `i'
}
drop incomebracket 
	
local balance_vars "elig_age elig_fpl female prior_enroll_2 zeropremium_eligible new_fpl new_elig"

mvreg `balance_vars' = arm2 arm3 arm4 if use_id==1 & letter_mailed < mdy(7,1,2019)

manova `balance_vars' = arm2 arm3 arm4 if use_id==1 & letter_mailed < mdy(7,1,2019)

mvtest means `balance_vars'  if use_id==1 & letter_mailed < mdy(7,1,2019), by(arm) heterogeneous 

* Generate matrix for F-statistics and P-values
mat demo_fstats = J(7,1,.)
mat demo_pvals = J(7,1,.)

global balance_vars elig_age elig_fpl female prior_enroll_2 zeropremium_eligible new_fpl new_elig

forv i=1/7 {
tokenize $balance_vars
reg ``i'' arm2 arm3 arm4 if use_id==1 & letter_mailed < mdy(7,1,2019), r
test arm2 == arm3 == arm4 ==0
mat demo_fstats[`i',1]=r(F)
mat demo_pvals[`i',1]=r(p)
}

mat fstats=demo_fstats
mat pvals=demo_pvals

mat results=fstats,pvals
mat rownames results = $balance_vars
mat colnames results = F-Stat P-Value

* Export matrices with results
mat2txt, saving("../output/main_results/tables/summary_fstats.tex") matrix(results) title("Full Sample") replace

* To get wald est of joint sig
local balance_vars "elig_age elig_fpl female prior_enroll_2 zeropremium_eligible new_fpl new_elig"

mvreg `balance_vars' = arm2 arm3 arm4 if use_id_90==1 & letter_mailed < mdy(7,1,2019)

manova `balance_vars' = arm2 arm3 arm4 if use_id_90==1 & letter_mailed < mdy(7,1,2019)

mvtest means `balance_vars'  if use_id_90==1 & letter_mailed < mdy(7,1,2019), by(arm) heterogeneous 
end

program main_results
use ../output/connector_regression_prep, clear

egen med_age = median(elig_age)
encode zipstring, gen(zip_factor)
gen over35 = (elig_age>35 & !mi(elig_age))
gen overmed = (elig_age > med_age & !mi(elig_age))
gen duesooner = (time_to_formdue < 30)

local folder "asof_Dec2020"

local basic_controls "female new_elig prior_enroll_2 spanish_letter i.batch i.WrapRegion_f time_elapsed new_fpl i.incomebracket i.incomebracket#c.new_fpl c.elig_age#i.age_cat"

local enhanced_controls "female new_elig prior_enroll_2 spanish_letter i.batch i.WrapRegion_f time_elapsed new_fpl i.incomebracket i.incomebracket#c.new_fpl c.elig_age#i.age_cat i.zip3digit"

di "No controls"
qui: reg enroll90 arm2 arm3 arm4 i.batch if use_id_90, r
est store a
parmest, saving("../output/main_results/`folder'/enroll90_noctrl_fullsample.dta", replace)

di "No Controls: Arm 4 - Arm 3"
lincom arm4 - arm3

* Arm 2 vs. Arm 3
test arm2 = arm3
mat nc_23_f = r(F)
mat nc_23_p = r(p)
sca nc_23_f = nc_23_f[1,1]
sca nc_23_p = nc_23_p[1,1]
estadd sca f_23 = nc_23_f
estadd sca p_23 = nc_23_p

* Arm 2 vs. Arm 4
test arm2 = arm4
mat nc_24_f = r(F)
mat nc_24_p = r(p)
sca nc_24_f = nc_24_f[1,1]
sca nc_24_p = nc_24_p[1,1]
estadd sca f_24 = nc_24_f
estadd sca p_24 = nc_24_p

* Arm 3 vs. Arm 4
test arm3 = arm4
mat nc_34_f = r(F)
mat nc_34_p = r(p)
sca nc_34_f = nc_34_f[1,1]
sca nc_34_p = nc_34_p[1,1]
estadd sca f_34 = nc_34_f
estadd sca p_34 = nc_34_p

* Arm 2 vs. Arm 3 vs. Arm 4 vs. 0
test arm2 = arm3 = arm4 = 0
mat nc_2340_f = r(F)
mat nc_2340_p = r(p)
sca nc_2340_f = nc_2340_f[1,1]
sca nc_2340_p = nc_2340_p[1,1]
estadd sca f_2340 = nc_2340_f
estadd sca p_2340 = nc_2340_p

tabstat enroll90 if arm1==1 & use_id_90==1, stats(mean) save
mat mean = r(StatTotal)
sca mean = mean[1,1]
estadd sca mean = mean
	
di "Basic controls" 
qui: reg enroll90 arm2 arm3 arm4 `basic_controls' if use_id_90==1, r
est store b
parmest, saving("../output/main_results/`folder'/enroll90_basicctrl_fullsample.dta", replace)

di "Basic Controls: Arm 4 - Arm 3"
lincom arm4 - arm3

* Arm 2 vs. Arm 3
test arm2 = arm3
mat bc_23_f = r(F)
mat bc_23_p = r(p)
sca bc_23_f = bc_23_f[1,1]
sca bc_23_p = bc_23_p[1,1]
estadd sca f_23 = bc_23_f
estadd sca p_23 = bc_23_p

* Arm 2 vs. Arm 4
test arm2 = arm4
mat bc_24_f = r(F)
mat bc_24_p = r(p)
sca bc_24_f = bc_24_f[1,1]
sca bc_24_p = bc_24_p[1,1]
estadd sca f_24 = bc_24_f
estadd sca p_24 = bc_24_p

* Arm 3 vs. Arm 4
test arm3 = arm4
mat bc_34_f = r(F)
mat bc_34_p = r(p)
sca bc_34_f = bc_34_f[1,1]
sca bc_34_p = bc_34_p[1,1]
estadd sca f_34 = bc_34_f
estadd sca p_34 = bc_34_p

* Arm 2 vs. Arm 3 vs. Arm 4 vs. 0
test arm2 = arm3 = arm4 = 0
mat bc_2340_f = r(F)
mat bc_2340_p = r(p)
sca bc_2340_f = bc_2340_f[1,1]
sca bc_2340_p = bc_2340_p[1,1]
estadd sca f_2340 = bc_2340_f
estadd sca p_2340 = bc_2340_p

tabstat enroll90 if arm1==1 & use_id_90==1, stats(mean) save
mat mean = r(StatTotal)
sca mean = mean[1,1]
estadd sca mean = mean

di "Enhanced controls" 
qui: reg enroll90 arm2 arm3 arm4 `enhanced_controls' if use_id_90==1, r
est store c
parmest, saving("../output/main_results/`folder'/enroll90_enhancedctrl_fullsample.dta", replace)

di "Enhanced Controls: Arm 4 - Arm 3"
lincom arm4 - arm3

* Arm 2 vs. Arm 3
test arm2 = arm3
mat ec_23_f = r(F)
mat ec_23_p = r(p)
sca ec_23_f = ec_23_f[1,1]
sca ec_23_p = ec_23_p[1,1]
estadd sca f_23 = ec_23_f
estadd sca p_23 = ec_23_p

* Arm 2 vs. Arm 4
test arm2 = arm4
mat ec_24_f = r(F)
mat ec_24_p = r(p)
sca ec_24_f = ec_24_f[1,1]
sca ec_24_p = ec_24_p[1,1]
estadd sca f_24 = ec_24_f
estadd sca p_24 = ec_24_p

* Arm 3 vs. Arm 4
test arm3 = arm4
mat ec_34_f = r(F)
mat ec_34_p = r(p)
sca ec_34_f = ec_34_f[1,1]
sca ec_34_p = ec_34_p[1,1]
estadd sca f_34 = ec_34_f
estadd sca p_34 = ec_34_p

* Arm 2 vs. Arm 3 vs. Arm 4 vs. 0
test arm2 = arm3 = arm4 = 0
mat ec_2340_f = r(F)
mat ec_2340_p = r(p)
sca ec_2340_f = ec_2340_f[1,1]
sca ec_2340_p = ec_2340_p[1,1]
estadd sca f_2340 = ec_2340_f
estadd sca p_2340 = ec_2340_p

tabstat enroll90 if arm1==1 & use_id_90==1, stats(mean) save
mat mean = r(StatTotal)
sca mean = mean[1,1]
estadd sca mean = mean

di "Lasso controls"
local lasso_controls "0.prior_enroll_2#c.time_elapsed 0.prior_enroll_2#0.overmed 1.incomebracket#0.wave1 1.female#1.overmed 1.incomebracket#1.duesooner 1.new_elig#3.incomebracket 0.prior_enroll_2#384.zip_factor 0.prior_enroll_2#11.zip3digit 3.incomebracket#1.overmed 1.new_elig#1.incomebracket 0.spanish_letter#11.zip3digit 0.female#10.WrapRegion_f 0.over35#386.zip_factor 11.WrapRegion_f#c.new_fpl 1.duesooner#24.zip3digit 1.prior_enroll_2#1.wave1 1.wave1#291.zip_factor 0.overmed#18.zip3digit 0.wave1#575.zip_factor 1.over35#310.zip_factor 25.zip3digit#c.new_fpl 0.new_elig#11.WrapRegion_f 5.incomebracket#1.wave1 0.duesooner#97.zip_factor 0.spanish_letter#381.zip_factor 402.zip_factor 0.spanish_letter#402.zip_factor 5.incomebracket#c.time_elapsed 0.new_elig#55.zip_factor 1.prior_enroll_2#4.incomebracket 1.female#416.zip_factor 1.female#3.incomebracket 12.WrapRegion_f#1.duesooner 1.prior_enroll_2#0.duesooner 0.overmed#386.zip_factor 1.incomebracket#16.zip3digit 1.new_elig#581.zip_factor 1.female#0.new_elig 1.overmed#271.zip_factor 1.female#310.zip_factor 3.incomebracket#343.zip_factor 0.new_elig#2.incomebracket 0.overmed#6.zip_factor 0.prior_enroll_2#18.zip3digit" 

qui: reg enroll90 arm2 arm3 arm4 i.batch `lasso_controls' if use_id_90==1, r
est store d
parmest, saving("../output/main_results/asof_Dec2020/enroll90_lassoctrl_fullsample.dta", replace)

di "LASSO Controls: Arm 4 - Arm 3"
lincom arm4 - arm3

* Arm 2 vs. Arm 3
test arm2 = arm3
mat lc_23_f = r(F)
mat lc_23_p = r(p)
sca lc_23_f = lc_23_f[1,1]
sca lc_23_p = lc_23_p[1,1]
estadd sca f_23 = lc_23_f
estadd sca p_23 = lc_23_p

* Arm 2 vs. Arm 4
test arm2 = arm4
mat lc_24_f = r(F)
mat lc_24_p = r(p)
sca lc_24_f = lc_24_f[1,1]
sca lc_24_p = lc_24_p[1,1]
estadd sca f_24 = lc_24_f
estadd sca p_24 = lc_24_p

* Arm 3 vs. Arm 4
test arm3 = arm4
mat lc_34_f = r(F)
mat lc_34_p = r(p)
sca lc_34_f = lc_34_f[1,1]
sca lc_34_p = lc_34_p[1,1]
estadd sca f_34 = lc_34_f
estadd sca p_34 = lc_34_p

* Arm 2 vs. Arm 3 vs. Arm 4 vs. 0
test arm2 = arm3 = arm4 = 0
mat lc_2340_f = r(F)
mat lc_2340_p = r(p)
sca lc_2340_f = lc_2340_f[1,1]
sca lc_2340_p = lc_2340_p[1,1]
estadd sca f_2340 = lc_2340_f
estadd sca p_2340 = lc_2340_p

tabstat enroll90 if arm1==1 & use_id_90==1, stats(mean) save
mat mean = r(StatTotal)
sca mean = mean[1,1]
estadd sca mean = mean

* Table
esttab a b c d, se star(+ 0.1 * 0.05 ** 0.01 *** 0.001) label keep(arm2 arm3 arm4 new_fpl female prior_enroll_2 spanish_letter _cons)
esttab a b c d using "../output/main_results/tables/`folder'/enroll90_fullsample.tex", ///
se ar2 lines replace b(%9.3fc) se(%9.3fc) star(+ 0.1 * 0.05 ** 0.01) label keep(arm2 arm3 arm4 new_fpl female prior_enroll_2 spanish_letter) scalars("mean Control Mean" "p_23 Arm 2 = Arm 3" "p_24 Arm 2 = Arm 4" "p_34 Arm 3 = Arm 4" "p_2340 Arm 2 = Arm 3 = Arm 4 = 0") ///	
title("Enroll by 90 Days (Full Sample)") ///
mtitle("No controls" "Basic controls" "Enhanced controls" "Post-LASSO controls") tex

end

* Appendix Tables
program mm_by_zp
use ../output/persistence_months, clear
egen first_12 = rowtotal(enr_1months_postletter-enr_13months_postletter)
egen first_18 = rowtotal(enr_1months_postletter-enr_19months_postletter)
egen first_24 = rowtotal(enr_1months_postletter-enr_25months_postletter)

tabstat first_12 first_18 first_24, stats(mean max)

foreach var in "first_12" "first_18" {
	if "`var'"=="first_12" local x = 12
	if "`var'"=="first_18" local x = 18
	
estimates clear
local folder "asof_Dec2020"	
* Full Sample
di "Outcome: `var', no controls, sample: `sample'" 
qui: reg `var' arm2 arm3 arm4 i.batch if use_id_90==1, r
est store a
parmest, saving("../output/main_results/`folder'/first_`x'_full.dta", replace)

* Arm 2 vs. Arm 3
test arm2 = arm3
mat full_23_f = r(F)
mat full_23_p = r(p)
sca full_23_f = full_23_f[1,1]
sca full_23_p = full_23_p[1,1]
estadd sca f_23 = full_23_f
estadd sca p_23 = full_23_p

* Arm 2 vs. Arm 4
test arm2 = arm4
mat full_24_f = r(F)
mat full_24_p = r(p)
sca full_24_f = full_24_f[1,1]
sca full_24_p = full_24_p[1,1]
estadd sca f_24 = full_24_f
estadd sca p_24 = full_24_p

* Arm 3 vs. Arm 4
test arm3 = arm4
mat full_34_f = r(F)
mat full_34_p = r(p)
sca full_34_f = full_34_f[1,1]
sca full_34_p = full_34_p[1,1]
estadd sca f_34 = full_34_f
estadd sca p_34 = full_34_p

* Arm 2 vs. Arm 3 vs. Arm 4 vs. 0
test arm2 = arm3 = arm4 = 0
mat full_2340_f = r(F)
mat full_2340_p = r(p)
sca full_2340_f = full_2340_f[1,1]
sca full_2340_p = full_2340_p[1,1]
estadd sca f_2340 = full_2340_f
estadd sca p_2340 = full_2340_p

tabstat `var' if arm1==1 & use_id_90==1, stats(mean) save
mat mean_full = r(StatTotal)
sca mean_full = mean_full[1,1]
estadd sca mean = mean_full

* ZP Elig.
qui: reg `var' arm2 arm3 arm4 i.batch if new_elig==1 & use_id_90==1, r
est store b
parmest, saving("../output/main_results/`folder'/first_`x'_zp1.dta", replace)

* Arm 2 vs. Arm 3
test arm2 = arm3
mat zp1_23_f = r(F)
mat zp1_23_p = r(p)
sca zp1_23_f = zp1_23_f[1,1]
sca zp1_23_p = zp1_23_p[1,1]
estadd sca f_23 = zp1_23_f
estadd sca p_23 = zp1_23_p

* Arm 2 vs. Arm 4
test arm2 = arm4
mat zp1_24_f = r(F)
mat zp1_24_p = r(p)
sca zp1_24_f = zp1_24_f[1,1]
sca zp1_24_p = zp1_24_p[1,1]
estadd sca f_24 = zp1_24_f
estadd sca p_24 = zp1_24_p

* Arm 3 vs. Arm 4
test arm3 = arm4
mat zp1_34_f = r(F)
mat zp1_34_p = r(p)
sca zp1_34_f = zp1_34_f[1,1]
sca zp1_34_p = zp1_34_p[1,1]
estadd sca f_34 = zp1_34_f
estadd sca p_34 = zp1_34_p

* Arm 2 vs. Arm 3 vs. Arm 4 vs. 0
test arm2 = arm3 = arm4 = 0
mat zp1_2340_f = r(F)
mat zp1_2340_p = r(p)
sca zp1_2340_f = zp1_2340_f[1,1]
sca zp1_2340_p = zp1_2340_p[1,1]
estadd sca f_2340 = zp1_2340_f
estadd sca p_2340 = zp1_2340_p

tabstat `var' if arm1==1 & new_elig==1 & use_id_90==1, stats(mean) save
mat mean_zp1 = r(StatTotal)
sca mean_zp1 = mean_zp1[1,1]
estadd sca mean = mean_zp1

*ZP Inelig.
qui: reg `var' arm2 arm3 arm4 i.batch if new_elig==0 & use_id_90==1, r
est store c
parmest, saving("../output/main_results/`folder'/first_`x'_zp0.dta", replace)

* Arm 2 vs. Arm 3
test arm2 = arm3
mat zp0_23_f = r(F)
mat zp0_23_p = r(p)
sca zp0_23_f = zp0_23_f[1,1]
sca zp0_23_p = zp0_23_p[1,1]
estadd sca f_23 = zp0_23_f
estadd sca p_23 = zp0_23_p

* Arm 2 vs. Arm 4
test arm2 = arm4
mat zp0_24_f = r(F)
mat zp0_24_p = r(p)
sca zp0_24_f = zp0_24_f[1,1]
sca zp0_24_p = zp0_24_p[1,1]
estadd sca f_24 = zp0_24_f
estadd sca p_24 = zp0_24_p

* Arm 3 vs. Arm 4
test arm3 = arm4
mat zp0_34_f = r(F)
mat zp0_34_p = r(p)
sca zp0_34_f = zp0_34_f[1,1]
sca zp0_34_p = zp0_34_p[1,1]
estadd sca f_34 = zp0_34_f
estadd sca p_34 = zp0_34_p

* Arm 2 vs. Arm 3 vs. Arm 4 vs. 0
test arm2 = arm3 = arm4 = 0
mat zp0_2340_f = r(F)
mat zp0_2340_p = r(p)
sca zp0_2340_f = zp0_2340_f[1,1]
sca zp0_2340_p = zp0_2340_p[1,1]
estadd sca f_2340 = zp0_2340_f
estadd sca p_2340 = zp0_2340_p

tabstat `var' if arm1==1 & new_elig==0 & use_id_90==1, stats(mean) save
mat mean_zp0 = r(StatTotal)
sca mean_zp0 = mean_zp0[1,1]
estadd sca mean = mean_zp0

* Table
di "Outcome: `lab', Sample: `sample'"
esttab a b c, se star(+ 0.1 * 0.05 ** 0.01 *** 0.001) label keep(arm2 arm3 arm4)
esttab a b c using "../output/main_results/tables/`folder'/first_`x'_zp.tex", ///
se lines replace b(%9.3fc) se(%9.3fc) star(+ 0.1 * 0.05 ** 0.01) label keep(arm2 arm3 arm4) scalars("mean Control Mean" "p_23 Arm 2 = Arm 3" "p_24 Arm 2 = Arm 4" "p_34 Arm 3 = Arm 4" "p_2340 Arm 2 = Arm 3 = Arm 4 = 0") ///	
title("`lab' (`sample')") ///
mtitle("Full Sample" "ZP Elig." "ZP Inelig.") tex
}
end

program het_effects
local folder "asof_Dec2020"
local population "fullsample overmed undermed zeroprem_elig zeroprem_inelig spanish english duesooner duelater"

local fullsample "use_id_90==1"
local zeroprem_elig "use_id_90==1 & new_elig==1"
local zeroprem_inelig "use_id_90==1 & new_elig==0"
local spanish "use_id_90==1 & spanish_letter==1"
local english "use_id_90==1 & spanish_letter==0"
local undermed "use_id_90==1 & elig_age <= med_age & elig_age!=."
local overmed "use_id_90==1 & elig_age > med_age & elig_age!=."
local duesooner "use_id_90==1 & time_to_formdue < 30" 
local duelater "use_id_90==1 & time_to_formdue > 30 & time_to_formdue!=."

local b = 0
foreach sample of local population {
	
use ../output/connector_regression_prep, clear
egen med_age = median(elig_age)
encode zipstring, gen(zip_factor)
	local b = `b'+1

di "Outcome: `var', no controls, sample: `sample'" 
qui: reg enroll90 arm2 arm3 arm4 i.batch if ``sample'', r
est store a`b'
parmest, saving("../output/main_results/`folder'/het_`sample'.dta", replace)

* Arm 2 vs. Arm 3
test arm2 = arm3
mat n`b'_23_f = r(F)
mat n`b'_23_p = r(p)
sca n`b'_23_f = n`b'_23_f[1,1]
sca n`b'_23_p = n`b'_23_p[1,1]
estadd sca f_23 = n`b'_23_f
estadd sca p_23 = n`b'_23_p

* Arm 2 vs. Arm 4
test arm2 = arm4
mat n`b'_24_f = r(F)
mat n`b'_24_p = r(p)
sca n`b'_24_f = n`b'_24_f[1,1]
sca n`b'_24_p = n`b'_24_p[1,1]
estadd sca f_24 = n`b'_24_f
estadd sca p_24 = n`b'_24_p

* Arm 3 vs. Arm 4
test arm3 = arm4
mat n`b'_34_f = r(F)
mat n`b'_34_p = r(p)
sca n`b'_34_f = n`b'_34_f[1,1]
sca n`b'_34_p = n`b'_34_p[1,1]
estadd sca f_34 = n`b'_34_f
estadd sca p_34 = n`b'_34_p

* Arm 2 vs. Arm 3 vs. Arm 4 vs. 0
test arm2 = arm3 = arm4 = 0
mat n`b'_2340_f = r(F)
mat n`b'_2340_p = r(p)
sca n`b'_2340_f = n`b'_2340_f[1,1]
sca n`b'_2340_p = n`b'_2340_p[1,1]
estadd sca f_2340 = n`b'_2340_f
estadd sca p_2340 = n`b'_2340_p

tabstat enroll90 if arm1==1 & ``sample'', stats(mean) save
mat mean_`b' = r(StatTotal)
sca mean_`b' = mean_`b'[1,1]
estadd sca mean = mean_`b'

}

esttab a1 a2 a3 a4 a5 a6 a7 a8 a9, se star(+ 0.1 * 0.05 ** 0.01 *** 0.001) label keep(arm2 arm3 arm4)
esttab a1 a2 a3 a4 a5 a6 a7 a8 a9 using "../output/main_results/tables/`folder'/het_effects.tex", ///
se lines replace b(%9.3fc) se(%9.3fc) star(+ 0.1 * 0.05 ** 0.01) label keep(arm2 arm3 arm4) scalars("mean Control Mean" "p_23 Arm 2 = Arm 3" "p_24 Arm 2 = Arm 4" "p_34 Arm 3 = Arm 4" "p_2340 Arm 2 = Arm 3 = Arm 4 = 0") ///	
mtitle("Overall" "Above Med. Age" "Below Med. Age" "ZP. Eligible" "ZP. Ineligible" "Spanish Letter" "English Letter" "Due Sooner" "Due Later") tex
end

program sum_stats_mail
import excel ${mail_comply}, firstrow clear
gen mail_sender = 1
duplicates drop
save ../output/mail_ids, replace

use ../output/jpal_sample, clear
keep if use_id==1 & letter_mailed < mdy(7,1,2019)

merge m:1 UniqueId using ../output/mail_ids, nogen keep(1 3)
replace mail_sender = 0 if mi(mail_sender)

* Separate by mail returner and enrollees in arm 4
gen group = .
replace group = 1 if mail_sender == 1 & enroll90==1 & arm==4
replace group = 2 if mail_sender == 1 & enroll90==0 & arm==4
replace group = 3 if mail_sender == 0 & enroll90==1 & arm==4
replace group = 4 if mail_sender == 0 & enroll90==0 & arm==4

* Summary statistics 
eststo clear
* Everyone in arm 4    
eststo: estpost summarize elig_age female prior_enroll_2 new_fpl new_elig if arm==4

* Returners vs. Not
sort mail_sender
* Did not return
eststo: estpost summarize ///
    elig_age female prior_enroll_2 new_fpl new_elig if arm==4 & mail_sender==0
* Did return
    eststo: estpost summarize ///
    elig_age female prior_enroll_2 new_fpl new_elig if arm==4 & mail_sender==1
    
* By mail returner and enroll by 90
sort group
by group: eststo: estpost summarize ///
    elig_age female prior_enroll_2 new_fpl new_elig if arm==4
    
esttab using "../output/main_results/tables/mail_summarytable.tex", replace ///
cells("mean sd")  label ///
title("Descriptive statistics for Arm 4") ///
mtitle("Full Arm 4" "Did Not Return - Total" "Did Return - Total" "Did Return - Enrolled by 90" "Did Return - Did Not Enroll by 90" "Did Not Return - Enrolled by 90" "Did Not Return - Did Not Enroll by 90") tex


* F-Stat and P-Value Matrices
mat arm1_fstats = J(5,1,.)
mat arm1_pvals = J(5,1,.)
mat arm3_fstats = J(5,1,.)
mat arm3_pvals = J(5,1,.)

global balance_vars elig_age female prior_enroll_2 new_fpl new_elig

gen group1 = (group==1)
gen group2 = (group==2)
gen group3 = (group==3)
gen group4 = (group==4)

forv i=1/5 {
tokenize $balance_vars
reg ``i'' arm1 arm4 if use_id==1 & letter_mailed < mdy(7,1,2019), r
test arm1==arm4==0
mat arm1_fstats[`i',1]=r(F)
mat arm1_pvals[`i',1]=r(p)
}
mat fstats_1=arm1_fstats
mat pvals_1=arm1_pvals

forv i=1/5 {
tokenize $balance_vars
reg ``i'' arm3 arm4 if use_id==1 & letter_mailed < mdy(7,1,2019), r
test arm3==arm4==0
mat arm3_fstats[`i',1]=r(F)
mat arm3_pvals[`i',1]=r(p)
}

mat fstats_3=arm3_fstats
mat pvals_3=arm3_pvals

mat results=fstats_1,pvals_1,fstats_3,pvals_3
mat rownames results = $balance_vars
mat colnames results = F-Stat-1 P-Value-1 F-Stat-3 P-Value-3

* Export matrices with results
mat2txt, saving("../output/main_results/tables/mail_returner_fstats.tex") matrix(results) title("Mail Returner Sample") replace

* To get wald est of joint sig
local balance_vars "elig_age female prior_enroll_2 new_fpl new_elig"

* Arm 1 Wald
mvreg `balance_vars' = arm1 arm4 if use_id_90==1 & letter_mailed < mdy(7,1,2019)

manova `balance_vars' = arm1 arm4 if use_id_90==1 & letter_mailed < mdy(7,1,2019)

mvtest means `balance_vars'  if use_id_90==1 & letter_mailed < mdy(7,1,2019), by(mail_sender) heterogeneous 
di "Arm 1"

* Arm 3 Wald
mvreg `balance_vars' = arm3 arm4 if use_id_90==1 & letter_mailed < mdy(7,1,2019)

manova `balance_vars' = arm3 arm4 if use_id_90==1 & letter_mailed < mdy(7,1,2019)

mvtest means `balance_vars'  if use_id_90==1 & letter_mailed < mdy(7,1,2019), by(mail_sender) heterogeneous 
di "Arm 3"
end

program rd_premium
foreach lim of numlist 200 250 {
use ../output/rd_setup_first_enroll, clear

* Bins vs 50% Banwidth
foreach i of numlist 1 2 5 {
	preserve
	collapse (mean) FamilyLevelBalancePremiumAmt fpl, by(bin_`i')
	
	if `i'==1 local x = "replace"
	if `i'!=1 local x = " "
	gen new_fpl = fpl-`lim'
	
* Banwidth 50
gen distance_`lim' = `lim'-fpl
gen above_cutoff = (fpl>=`lim')
gen interact = new_fpl*above_cutoff
local lb = `lim'-50
local ub = `lim'+50

label var new_fpl "FPL (Centered)"
label var above_cutoff "Above Cutoff"
label var interact "Above Cutoff x FPL (Centered)"

reg FamilyLevelBalancePremiumAmt above_cutoff new_fpl interact if abs(distance_`lim')<=50, r
outreg2 using "../output/main_results/tables/rd_netprem_`lim'fpl.xls", `x' addtext(Bins, `i'%, Banwidth, 50% FPL) label ctitle("Ave. Net-of-Subsidy Premium") addnote(Note: Cutoff = `lim'% FPL, Sample: First Enrollment Period per UniqueId (2016-2019))

drop above_cutoff interact
restore
}

* 5% Bins vs. Banwidths
preserve
collapse (mean) FamilyLevelBalancePremiumAmt fpl, by(bin_5)
	
gen new_fpl = fpl-`lim'	
gen distance_`lim' = `lim'-fpl
gen above_cutoff = (fpl>=`lim')
gen interact = new_fpl*above_cutoff

label var new_fpl "FPL (Centered)"
label var above_cutoff "Above Cutoff"
label var interact "Above Cutoff x FPL (Centered)"

foreach bw of numlist 40 30 20 {
reg FamilyLevelBalancePremiumAmt above_cutoff new_fpl interact if abs(distance_`lim')<=`bw', r
outreg2 using "../output/main_results/tables/rd_netprem_`lim'fpl.xls", addtext(Bins, 5%, Banwidth, `bw'% FPL) label ctitle("Ave. Net-of-Subsidy Premium")
}
restore

* 5% Bin vs 50% Banwidth - FPL-Squared
preserve
collapse (mean) FamilyLevelBalancePremiumAmt fpl, by(bin_5)

gen new_fpl = fpl-`lim'
gen distance_`lim' = `lim'-fpl
gen above_cutoff = (fpl>=`lim')
gen interact = new_fpl*above_cutoff
gen new_fpl_2 = new_fpl*new_fpl
gen interact_2 = new_fpl_2*above_cutoff

label var new_fpl "FPL (Centered)"
label var above_cutoff "Above Cutoff"
label var interact "Above Cutoff x FPL (Centered)"
label var new_fpl_2 "FPL-Squared (Centered)"
label var interact_2 "Above Cutoff x FPL-Squared (Centered)"

reg FamilyLevelBalancePremiumAmt above_cutoff new_fpl interact new_fpl_2 interact_2 if abs(distance_`lim')<=50, r
outreg2 using "../output/main_results/tables/rd_netprem_`lim'fpl.xls", addtext(Bins, 5%, Banwidth, 50% FPL) label ctitle("Ave. Net-of-Subsidy Premium")
restore
}
end

program rd_log
foreach lim of numlist 200 250 {
use ../output/rd_setup_first_enroll, clear

* 200 FPL
foreach i of numlist 1 2 5 {
	preserve
	collapse (sum) n (mean) fpl bin_2 bin_5, by(bin_1)
	collapse (mean) n fpl, by(bin_`i')
	
	if `i'==1 local x = "replace"
	if `i'!=1 local x = " "
	gen new_fpl = fpl-`lim'
	
* Banwidth 50
gen distance_`lim' = `lim'-fpl
gen above_cutoff = (fpl>=`lim')
gen interact = new_fpl*above_cutoff

label var new_fpl "FPL (Centered)"
label var above_cutoff "Above Cutoff"
label var interact "Above Cutoff x FPL (Centered)"

gen ln_enroll = log(n+1)

reg ln_enroll above_cutoff new_fpl interact if abs(distance_`lim')<=50, r
local num = _b[above_cutoff]
local den = _b[_cons]
local elasticity = `num'/`den'
local elasticity = round(`elasticity',0.01)

outreg2 using "../output/main_results/tables/rd_log_`lim'fpl.xls", `x' addtext(Elasticity, `elasticity', Bins, `i'%, Banwidth, 50% FPL) label ctitle("Log New Enroll") addnote(Note: Cutoff = `lim'% FPL, Sample: First Enrollment Period per UniqueId (2016-2019))

drop above_cutoff interact
restore
}

preserve
collapse (sum) n (mean) fpl bin_2 bin_5, by(bin_1)
collapse (mean) n fpl, by(bin_5)
	
gen new_fpl = fpl-`lim'	
gen distance_`lim' = `lim'-fpl
gen above_cutoff = (fpl>=`lim')
gen interact = new_fpl*above_cutoff

label var new_fpl "FPL (Centered)"
label var above_cutoff "Above Cutoff"
label var interact "Above Cutoff x FPL (Centered)"

gen ln_enroll = log(n+1)

foreach bw of numlist 40 30 20 {
reg ln_enroll above_cutoff new_fpl interact if abs(distance_`lim')<=`bw', r
local num = _b[above_cutoff]
local den = _b[_cons]
local elasticity = `num'/`den'
local elasticity = round(`elasticity',0.01)

outreg2 using "../output/main_results/tables/rd_log_`lim'fpl.xls", addtext(Elasticity, `elasticity', Bins, 5%, Banwidth, `bw'% FPL) label ctitle("Log New Enroll")

}
restore

preserve
collapse (sum) n (mean) fpl bin_2 bin_5, by(bin_1)
collapse (mean) n fpl, by(bin_5)

gen new_fpl = fpl-`lim'
gen distance_`lim' = `lim'-fpl
gen above_cutoff = (fpl>=`lim')
gen interact = new_fpl*above_cutoff
gen new_fpl_2 = new_fpl*new_fpl
gen interact_2 = new_fpl_2*above_cutoff

label var new_fpl "FPL (Centered)"
label var above_cutoff "Above Cutoff"
label var interact "Above Cutoff x FPL (Centered)"
label var new_fpl_2 "FPL-Squared (Centered)"
label var interact_2 "Above Cutoff x FPL-Squared (Centered)"

gen ln_enroll = log(n+1)

reg ln_enroll above_cutoff new_fpl interact new_fpl_2 interact_2 if abs(distance_`lim')<=50, r
local num = _b[above_cutoff]
local den = _b[_cons]
local elasticity = `num'/`den'
local elasticity = round(`elasticity',0.01)

outreg2 using "../output/main_results/tables/rd_log_`lim'fpl.xls", addtext(Elasticity, `elasticity', Bins, 5%, Banwidth, 50% FPL) label ctitle("Log New Enroll")
restore
}
end

* Execute
*main
